//
//  Mat.h
//  Cartoonization
//
//  Created by Adrien MATRICON on 05/03/13.
//  Copyright (c) 2013 Adrien MATRICON. All rights reserved.
//

#ifndef Cartoonization_Mat_h
#define Cartoonization_Mat_h

#include <vector>
#include <iostream>

struct Mat
{
	std::vector<double> vals;
	int n, m;
	
	Mat(){}
	Mat(int a_n, int a_m): vals(a_n*a_m, 0), n(a_n), m(a_m) {}
	Mat(const Mat &mat): vals(mat.n * mat.m), n(mat.n), m(mat.m)
	{
		for(int i=0; i<n; i++)
		{
			for(int j=0; j<m; j++)
			{
				vals[i*m+j] = mat.vals[i*m+j];
			}
		}
	}
	
	Mat(const Vec v): vals(3), n(3), m(1)
	{
		vals[0] = v.x;
		vals[1] = v.y;
		vals[2] = v.z;
	}
	
	inline double operator()(int i, int j) const
	{
		return vals[i*m+j];
	};
	
	inline double& operator()(int i, int j)
	{
		return vals[i*m+j];
	};
	
	inline Mat transpose()
	{
		Mat result(m,n);
		
		for(int i=0; i<n; i++)
		{
			for(int j=0; j<m; j++)
			{
				result(j,i) = vals[i*m+j];
			}
		}
		
		return(result);
	}
	
	inline Mat& operator+= (const Mat &mat)
	{
		if(n != mat.n || m!= mat.m)
		{
			std::cerr<<"Matrix sizes do not match in addition method"<<std::endl;
			exit(1);
		}
		
		for(int i=0; i<n*m; i++)
		{
			vals[i] += mat.vals[i];
		}
		return(*this);
	};
	
	inline Mat& operator-= (const Mat &mat)
	{
		if(n != mat.n || m!= mat.m)
		{
			std::cerr<<"Matrix sizes do not match in soustraction method"<<std::endl;
			exit(1);
		}
		
		for(int i=0; i<n*m; i++)
		{
			vals[i] -= mat.vals[i];
		}
		return(*this);
	};
	
	inline Mat& operator*= (double s)
	{
		for(int i=0; i<n*m; i++)
		{
			vals[i] *= s;
		}
		return(*this);
	};
	
	inline Mat& operator/= (double s)
	{
		for(int i=0; i<n*m; i++)
		{
			vals[i] /= s;
		}
		return(*this);
	};
	
	//Returns the solution of A*x=B
	static Mat gauss(Mat A, Mat B)
	{
		if(A.n != B.n)
		{
			std::cerr<<"Matrix sizes do not match in gauss resolution"<<std::endl;
			exit(1);
		}
		
		int n = A.n;
		
		for(int i=0; i< n; i++)
		{
			if(A(i,i) == 0)
			{
				int counter = i+1;
				
				while (A(counter,i) == 0 && counter < n)
				{
					counter++;
				}
				
				if(counter == n)
				{
					std::cerr<<"This is not a Cramer system"<<std::endl;
					exit(1);
				}
				
				double tmp;
				
				for(int j=0; j<n; j++)
				{
					tmp = A(i,j);
					A(i,j) = A(counter,j);
					A(counter,j) = tmp;
				}
			
				tmp = B(i,0);
				B(i,0) = B(counter,0);
				B(counter,0) = tmp;
			}

			double tmp = A(i,i);

			for(int j=0; j<n; j++)
			{
				A(i,j) /= tmp;
			}

			B(i,0) /= tmp;

			for(int k = i+1; k<n; k++)
			{
				tmp = A(k,i);
				
				if(tmp != 0)
				{
					for(int j=0; j<n; j++)
					{
						A(k,j) -= tmp*A(i,j);
					}
					
					B(k,0) -= tmp*B(i,0);
				}
			}
		}

		Mat result(n,1);

		for(int i=n-1; i >= 0; i--)
		{
			result(i,0) = B(i,0);
			
			for(int j=i+1; j < n; j++)
			{
				result(i,0) -= A(i,j)*result(j,0);
			}
			
			result(i,0) /= A(i,i);
		}

		return(result);
	}
};


inline Mat operator* (const Mat &m1, const Mat &m2)
{
	if(m1.m != m2.n)
	{
		std::cerr<<"Matrix sizes do not match in product"<<std::endl;
		exit(1);
	}
	
	Mat result (m1.n, m2.m);
	
	for(int i=0; i<m1.n; i++)
	{
		for(int j=0; j<m2.m; j++)
		{
			for(int k=0; k<m1.m; k++)
			{
				result(i,j) += m1(i,k) * m2(k,j);
			}
		}
	}
	
	return(result);
};

inline Mat operator+ (const Mat &m1, const Mat &m2)
{
	if(m1.n != m2.n || m1.m != m2.m)
	{
		std::cerr<<"Matrix sizes do not match in addition"<<std::endl;
		exit(1);
	}
	
	Mat result (m1.n, m1.m);
	
	for(int i=0; i<m1.n; i++)
	{
		for(int j=0; j<m2.m; j++)
		{
			result(i,j) = m1(i,j) + m2(i,j);
		}
	}
	
	return(result);
}

inline Mat operator- (const Mat &m1, const Mat &m2)
{
	if(m1.n != m2.n || m1.m != m2.m)
	{
		std::cerr<<"Matrix sizes do not match in soustraction"<<std::endl;
		exit(1);
	}
	
	Mat result (m1.n, m1.m);
	
	for(int i=0; i<m1.n; i++)
	{
		for(int j=0; j<m2.m; j++)
		{
			result(i,j) = m1(i,j) - m2(i,j);
		}
	}
	
	return(result);
}

inline Mat operator* (const Mat &mat, const double &s)
{
	Mat result (mat.n, mat.m);
	
	for(int i=0; i<mat.n; i++)
	{
		for(int j=0; j<mat.m; j++)
		{
			result(i,j) = mat(i,j) * s;
		}
	}
	
	return(result);
};

inline Mat operator* (const double &s, const Mat &mat)
{
	Mat result (mat.n, mat.m);
	
	for(int i=0; i<mat.n; i++)
	{
		for(int j=0; j<mat.m; j++)
		{
			result(i,j) = mat(i,j) * s;
		}
	}
	
	return(result);
};

inline Mat operator/ (const Mat &mat, const double &s)
{
	Mat result (mat.n, mat.m);
	
	for(int i=0; i<mat.n; i++)
	{
		for(int j=0; j<mat.m; j++)
		{
			result(i,j) = mat(i,j) / s;
		}
	}
	
	return(result);
};

#endif
